SCIENTIFIC 

REPORTS 



s 




OPEN 



SUBJECT AREAS: 

ELECTRONIC STRUCTURE 

ELECTRONIC PROPERTIES AND 
MATERIALS 

APPLIED PHYSICS 

ELECTRICAL AND ELECTRONIC 
ENGINEERING 



The transport properties of oxygen 
vacancy-related polaron-like bound state 
in HfO, 



X 



Zhongrui Wang\ HongYu Yu^ & Haibin Su^ 



^ Nanyang Technological University, Singapore 639798, ■^South University of Science and Technology of China, China 5 1 8055. 



Received 
2 July 2013 

Accepted 
1 November 201 3 

Published 
9 December 2013 



Correspondence and 
requests for materials 
should be addressed to 
H.Y.Y. (yu.hy@sustc. 
edu.cn) or H.B.S. 
(hbsu@ntu.edu. sg) 



The oxygen vacancy-related polaron-like bound state migration in HfOx accounting for the observed 
transport properties in the high resistance state of resistive switching is investigated by the density 
functional theory with hybrid functional. The barrier of hopping among the threefold oxygen vacancies is 
strongly dependent on the direction of motion. Especially, the lowest barrier along the <001> direction is 
90 meV, in agreement with the experimental value measured from 135 K to room temperature. This 
hopping mainly invokes the z-directional motion of hafnium and threefold oxygen atoms in the vicinity of 
the oxygen vacancy resulted from the synergized combination of coupled phonon modes. In the presence of 
(111) surface, the lowest barrier of hopping between the surface oxygen vacancies is 360 meV along the 
< 101 > direction, where the significant surface perpendicular motion of hafnium and twofold oxygen atoms 
surrounding the oxygen vacancy is identified to facilitate this type of polaron-like bound state migration. 
Thus, the migration on the surfaces could be more important at the high temperature. 

Resistive switching in thin film devices attracts great interests for potential applications in future non -volatile 
memories The monoclinic phase of HfOx (m-HfOx) is a material of particular importance to resistive 
switching memories, which exhibits reversible metal-insulator transitions in operations''^. Oxygen vacan- 
cies have been widely speculated to be responsible for the valence change type bipolar resistive switching in 
oxides^'^. Spectroscopic ellipsometer measurements have revealed that oxygen vacancies are dominant native 
defects in m-HfOx^. The observed switching process in m-HfOx is likely due to the percolation of oxygen 
vacancies^. In the high resistance state, the electronic state associated with doped oxygen vacancies tends to 
localize within the mobility gap, which possesses experimental significances of polaron-like transport as revealed 
by the field and temperature dependent dielectric loss measurements^. And the observed thermal activation 
energy transition at half of the Debye temperature (Ea ~ 50 meV at T > 0/2) is consistent with the prediction 
based on the kinetics of polaronic motion by Schnakenberg^"^^. 

The presence of polarons in transition metal oxides are facilitated by the significant electron-phonon 
interaction^ The details of polaron transport are reviewed by Holstein^^'^^, Marcus^^, Emin^^, Friedman^^, 
Alexandrov^^, Austin and Mott^°, Ross et al.^^ The intrinsic polaron state has been proposed in m-Hf02 without 
invoking oxygen vacancies^°'^\ For oxygen deficient HfOx, the trapping of extra electrons by defect states 
promotes electron -lattice coupling, which leads to the formation of polaron-like bound states in m-HfOx^^. 
Density functional calculations of polaron-like bound state trapping on oxygen vacancies of m-HfOx were 
reported with hybrid functionals approaches^'"^^ Recently, continuous calculation by using density functional 
theory of the total energy along the entire migration pathway has been reported of accurate predictions of 
activation energy^^'^^"^^. Despite intensive studies of polaron-like carrier trapping in m-HfOx, the detailed migra- 
tion pathway of oxygen vacancy- related polaron-like bound states in m-HfOx in either bulk lattice or presence of 
surfaces is yet to be clarified, which is of great importance for interpreting the observed transport behaviours in m- 
HfOx based resistive switching memory. Besides, the influence of lattice on the polaron-like bound state hopping 
is unclear in m-HfOx- Here we study the high resistance state transport in oxygen deficient m-HfOx resistive 
memory leveraged on the previous studies on charge migration^^'^^'^^'^^. Computational evidence is presented for 
the motion of oxygen vacancy-related polaron-like bound state in m-HfOx- Using hybrid functional based 
conjugate gradient refined linear and quadratic synchronous transit schemes'\ we characterize the anisotropic 
hopping barriers for polaron-like bound states jumping between equivalent threefold oxygen vacancies along 
lattice vectors. The value of the lowest barrier is in consistency with experimental measurement from 135 K to 
room temperature. Projected phonon density of states has been calculated, indicating large vibration with 
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Figure 1 | (a) The optimized 1X1X2 lattice with two threefold oxygen 
vacancies. Hafnium, oxygen atom, and vacancy are represented by blue, 
red and grey balls, respectively. An isosurface (Isovalue = 0.05 A"^) shows 
the spin distribution. The black dashed line and labels indicate inter- 
atomic distances in angstroms between the vacancy-neighbouring 
hafnium atoms. The spin mainly localizes on the lower oxygen vacancy site, 
(b) (100) and (001) slice view of the electron localization function 
calculated from valence electrons of the configuration in (a). The arrows 
point to the original positions of the oxygen vacancies. An attractor is 
located at the site of the lower vacancy while the value of electron 
localization function of the upper site is close to zero, (c) Atom projected 
density of states for the upper unit cell (non-polar site) and lower unit cell 
(polar site) in (a). Fermi level is set to be zero eV. The polaron-like bound 
state is mainly composed of Hf 5d atomic orbital. 

vacancy-neighbouring threefold oxygen and hafnium atoms. In the 
presence of (111) surface, the kinetic barrier of the symmetrical 
hopping is significantly higher than that of the hopping in the bulk. 
Surface phonon modes, localizing on vacancy-neighbouring haf- 
nium and twofold oxygen atoms, are identified to facilitate this hop- 
ping. The experimentally observed barrier in the high temperature 
range from 320 K to 400 K is clearly about two to three times larger 
than that in low temperature. Thus, this work suggests that surface 
polaron migration may be important to the observed transport at 
high temperature. 



Table 1 Calculated spin 


distribution over vacancy neighboring 


hafnium atoms 




Lower vacancy site (H) 


Upper vacancy site (H) 


0.11 


0.01 


0.16 


0.01 


0.22 


0 



Results 

The imperfect treatment of self-interaction in density functional 
theory tends to overly delocalize d electrons and yield an under- 
estimation of the bandgap^^. Recently functionals such as sX, PBEO 
and B3LYP have been used in the study of oxygen vacancies in 
HfOx^^"^^. By including a fraction of exact exchange, the B3LYP 
functional in our calculation, which correctly predicts stoichiometric 
monoclinic HfOs (a = 5.00 A, b = 5.06 A, c = 5.18 A, |3 = 100.1°) to 
be an insulator with a gap of 6.2 eV in good agreement with previous 
reports^'^^"^^. Introducing oxygen vacancies with positive charge cre- 
ates t2g molecular orbitals with a symmetry to hafnium-vacancy 
axis^^. Here, each supercell contains two m-Hf02 unit cells with 
oxygen vacancies at equivalent positions. Although the threefold 
oxygen vacancy has strong tendency of retaining +2 charge state^^"^^, 
the total charge per supercell is set to be +3 to model the effect of 
electron reservoirs existed during the operation of RRAM devices. 

Structure- wise, with 12.5% concentration of fourfold oxygen 
vacancies, the lattice is transformed into orthorhombic phase (a = 
10.58 A, b = 4.20 A, c = 4.27 A). Subsequent transition to tetragonal 
phase is observed with increasing non-stoichiometry, consistent with 
the report of Xue et al.^^ On the other hand, threefold oxygen vacan- 
cies only transform the structure into triclinic phase. The structural 
transformation resulting from the increase of the oxygen vacancy 
content is also observed in other transition metal oxides including 
cuprates^^ The localization feature of the polaron-like bound state 
needs judicious correction of self-interaction error. Localization of 
the polaron-like bound state is studied with the exchange limit on 
both types of oxygen vacancies. For fourfold oxygen vacancies, the 
charge density computed by exact exchange is localized on one oxy- 
gen vacancy while the charge is distributed over both vacancies when 
computed with B3LYP functional, which suggests that the local- 
ization in fourfold oxygen vacancies requires much stringent treat- 
ment of the self-interaction error. However, the charge is localized on 
threefold oxygen vacancies when computed with B3LYP, which is 
consistent with the result by exact exchange. Moreover, threefold 
oxygen vacancies are more energetically favourable. Thus, the main 
study focuses on threefold oxygen vacancies. 

The oxygen vacancy in m-HfOx induces changes of the local lattice 
structure. In the 12.5% concentration case, after B3LYP optimiza- 
tion, phase transformation takes place in the 1 X 1 X 2 supercell, 
which transforms the lattice to triclinic phase, (a = 4.71 A, b = 4.95 
A, c = 9.60 A, a = 90.9°, |3 = 95.2°, y = 89.7°) The localization of the 
polaron-like bound state shown in Figure 1(a) is accompanied by 
displacement of three neighbouring hafnium atoms by —0.2 A 
towards each other, on average, about 10% of the O-Hf distance. 
In this case, about 50% of the electron spin density is locaUzed pre- 
dominantly on three hafnium ions sharing the missing threefold 
oxygen atoms, similar to that of the polaron in m-Hf02^°. The spin 
distribution, listed in Table 1, is asymmetrical among the hafnium 
ions. Electron localization function has been calculated from the 
valence electrons as an approximation to the all-orbital electron 
localization function. (Note that electron localization function values 
from CASTEP do not span the full range from 0 to 1.) Irreducible /- 
localization domain, as shown in Figure 1(b), could be defined by 
comprising the single attractor located at the position of the lower 
oxygen vacancy. As a contrast, the calculated electron localization 
function in the vicinity of the higher oxygen vacancy is close to zero. 
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Table 2 | Calculated lowest polaron-like bound state migration 
barriers between threefold oxygen vacancies in the equivalent sites 
in the bulk HfOx along < 1 00>, <01 0> and <001 > directions 



Direction Hopping Barrier (meV) Distance (A) 



<100> 661.19 4.62 

<010> 2814.79 4.76 

<001> 90.47 4.80 



which confirms the asymmetrical spin locaUzation. Figure 1(b) 
shows the atomic projected density of states. The polaron-Uke bound 
state occupies the spin-up state at the Fermi level, and is mainly 
composed of hafnium 5d orbitals on which the electron is localized. 
The localization splits the original degenerate bandgap states 
induced by oxygen vacancies by pushing up the un-occupied spin- 
down state 2 eV above the Fermi level. Weak coupling between 
vacancies leads to the small band dispersion (—0.26 eV) of the 
polaron-like bound state, implying electron localization in the low 
density limit^^. 

The polaron-like bound state migration is strongly dependent on 
the direction of motion in bulk HfOx- Polaron-like bound states 
migrate along three different directions (i.e. <100>, <010>, and 
<001>) between two equivalent neighbouring vacancies. The hop- 
ping direction dependent barriers unravel the anisotropic character- 
istics in transport process as presented in Table 2. Austin and Mott 
have proposed that the relationship between polaron hopping activa- 

eW 1 l\fl l\ , 

tion energy and distance is Wh = — , where e 

4 V^oo Kj\rp Rj' 

is the electron charge, kJk is the high-frequency/static dielectric 
constant, Vp is the polaron radius, and R is the distance between 
hopping centres^". However, the application of this formula in aniso- 
tropic transport needs more caution. In Figure 2(a), the transition 
state of <100> directional migration is accompanied with remark- 
able lattice distortions. Appreciable changes of bond length are 
observed, as shown in Figure 2(a). Compared to the charge in the 
initial localized state, the charge in the same region is reduced when 
in the transition state. This leads to the increment of the inter- atomic 
separation between hafnium atoms neighbouring the lower oxygen 



vacancy where the charge is originally trapped. So, the chemical 
bonding between hafnium and oxygen has covalent component 
besides the commonly discussed ionic nature. Similar analysis is 
performed for the hopping along <010> and <001> directions. 
Hopping along <010> direction experiences the highest barrier as 
shown in Figure 2(b), where significant changes in bond length are 
identified between threefold oxygen and hafnium atoms. When hop- 
ping is along <001> direction, the changes in bond lengths are 
much smaller mainly involving atoms near the vacancies as shown 
in Figure 2(c). The computed lowest activation energy is 90 meV 
along the <001 > direction, which is in agreement with the measured 
thermal activation energy of 87 meV observed by Compagnoni et al. 
in the temperature range from 135 K to room temperature^^. Ramo et 
al. have reported that the electron polaron migration barrier in per- 
fect Hf02 is 50 meV by linear extrapolation approach^°. Compared 
to hopping of polaron, polaron-like bound state migration is facing a 
relatively larger barrier. 

We have calculated the projected phonon density of states, 
shown in Figure 3 (a), of the configuration with lowest hopping 
barrier along <001> direction. Projecting the displacement vec- 
tor, the coordinate difference between transition state and reactant 
factored by the square root of atomic mass, onto the phonon 
eigenvectors at T point, the eigenvalues of the modes with top 
ten projections span the energy range from 8 to 45 meV. By 
inspecting the eigenvectors, the vibration of the mode with largest 
projection is polarized perpendicular to z-axis, mainly concerning 
the motion of hafnium atoms and threefold oxygen atoms. Higher 
order effect of the phonon eigenvectors are also estimated by 
varying the amplitude of the displacements, where the modes 
could be grouped based on the inter-mode anharmonic coupling. 
The transition state displacement is mainly a superposition of two 
groups of phonon modes, one concerning the modes localized on 
hafnium atoms with energy from 8 to 14 meV while the other 
consisting of motion of threefold oxygen atoms parallel to z-axis 
(i.e. 03 and 04) with energy from 24 to 45 meV, as illustrated in 
the inset of Figure 3(a) and Figure 3(b). The synergized combina- 
tion of strongly coupled phonon modes play important roles in 
the migration may closely relate to the migration of polaron along 
the <001> direction. 



(a) \/ 




1^ 



5^ fi) 





X ^ Y ^-xl^ 



Figure 2 | Calculated transition states of (a) <100> directional hopping in a 2 XIX llattice with two threefold oxygen vacancies (b) <010> directional 
hopping in a 1 X 2 X 1 lattice with two threefold oxygen vacancies (c) <001> directional hopping in a 1 X 1 X 2 lattice with two threefold oxygen 
vacancies. The polaron-like bound states are originally localized on lower oxygen vacancies. The black dashed lines and labels indicate changes of inter- 
atomic distances in angstroms between the vacancy-neighbouring hafnium atoms. Green labels are major bond length changes given in angstroms. 
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Density of Phonon States (A. U.) 

Figure 3 | (a) Calculated atom projected phonon density of states of the 1X1X2 HfOx lattice in Figure 1. The lower inset shows the zoomed view from 
8 meV to 14 meV while the upper inset shows the zoomed view from 23 to 43 meV, which correspond to the contributions of the labelled threefold 
oxygen atoms in (b). (b) Illustration of combined phonon packets. Arrows are the superposition of the two dominant phonon groups of coupled modes 
where the vibration is mainly with hafnium and threefold oxygen atoms. 



McKenne and Shluger have found that grain boundaries act as 
sinks for oxygen vacancies so percolation paths are preferentially 
located along grain boundaries'^. Thus it is interesting to investigate 
the polaron-hke bound state migration pathways in the vicinity of the 
boundaries. In this work, we focus on the migration dynamics near 
the (111) surface which possesses the lowest surface energy (1.105 J/ 
m^) among various low index surfaces^. For simplicity, symmetrical 
hopping between equivalent sites on (111) surfaces are investigated 
along either <110> or <101> direction with sixteen different 
paths. By using the slab model with two layers of HfOx structure (a 
= 6.27 A, b = 13.76 A), the lowest activation barrier is found to be 
360 meV along the < 101 > direction over 6.24 A as given in Table 3. 
Experimentally, the trap depth in HfOx, which presumes a Poole- 
Frenkel emission model, was found to fall in the range from 210 meV 
to 350 meV in the elevated temperature region from 300 K to 400 
j^39,4o Poole-Frenkel behaviours are likely to have an origin of 
small polar on hopping^'^^ it is very likely that polaron-like bound 
state hopping near surface dominates over the bulk hopping in the 
high temperature region due to oxygen vacancy segregation near 
surfaces and the subsequent suppression of inter-grain coherent 
linkage. 

Figure 4 (a) is the projected phonon density of states. Projecting 
the displacement vector onto the phonon eigenvectors at F-point, 
the modes with largest projections are below 48 meV, especially in 
the range from 7 to 20 meV. The dominant contribution to branches 
with large projections comes from the vibration of hafnium and 
twofold oxygen atoms neighbouring oxygen vacancies, as shown in 
Figure 4(a). Analysis of the phonon eigenvector with largest projec- 
tion reveals the vibration is mainly polarized perpendicular to the 
surface. It is also noted that relatively large projection comes from 
branches with energy 45 meV to 48 meV, passing through the for- 
bidden regions between the bulk continuums at non-zero (^-points in 
Figure 4(b), which reveals the importance of surface modes to 
polaron hopping along the surface^^. These branches are also noticed 
to be almost dispersionless along the symmetry directions due to the 



very large mass difference between hafnium and oxygen atoms. We 
have also computed the phonon structure of the (111) surface with- 
out oxygen vacancies and strong resemblance is found in dispersion. 
An exception, however, is that the branches appearing in the forbid- 
den region between the acoustic and optical bulk phonon modes are 
drastically changed due to oxygen vacancies induced surface relaxa- 
tion. Like the bulk case, anharmonic coupling between modes is 
estimated by varying the amplitude used in the finite displacement 
method. The displacement of the transition states is mainly com- 
posed of the group of coupled modes in the energy range from 7 to 
13 meV, as shown in Figure 4(c), localized on the hafnium and 
twofold oxygen atoms in the vicinity of the oxygen vacancies. 

Discussion 

The oxygen vacancy- related polaron-like bound state is clearly 
localized in the high resistance state of m-HfOx investigated by 
the density functional theory with the hybrid functional including 
proper portion of non-local exact exchange. Anisotropic migra- 
tion barriers are found for polaron-like bound state hopping 
between the equivalent threefold oxygen vacancies. The lowest 
barrier of threefold oxygen vacancy- related polaron-like bound 
state along the <001> direction is 90 meV, in agreement with 
the experimental value. The z-directional motion of the neigh- 
bouring hafnium and threefold oxygen atoms are identified to 
play important roles in the hopping of polaron-like bound state. 
In the presence of (111) surface, lowest migration hopping path is 



Table 3 | Calculated lowest polaron-like bound state migration 
barriers between threefold oxygen vacancies in the equivalent sites 
on the (111) surface of HfOx along <1 10> and <101 > directions 

Direction Hopping Barrier (meV) Distance (A) 

<110> 519.70 6.88 

<101> 360.93 6.24 
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Figure 4 | (a) Calculated atom projected phonon density of states and (b) phonon dispersion of (111) HfOx 2-layer vacuum slab model with lowest 
symmetrical migration energy along <101> direction. The projected bulk phonon dispersion of HfOx is shown by the hatched region. The branches 
within energy range from 45 to 48 meV pass through the forbidden regions between the bulk continuums at non-zero ^-points, representing surface 
modes, (c) Illustration of the dominant phonon packet. Arrows are the superposition of coupled modes where the vibration is mainly with the vacancy- 
neighbouring hafnium and twofold oxygen atoms. 



found to be along <101> direction with displacements mainly 
invoking the surface perpendicular motion of phonon modes loca- 
lized on the vacancy-neighbouring hafnium and twofold oxygen 
atoms. The present work based on the oxygen vacancy- related 
polaron-like bound state sheds hght on the microscopic nature 
of the transport process in the high resistance state of m-HfOx 
resistive memory devices. 

Method 

Here we use the plane wave pseudo-potential code CASTEP"^^, with norm-conserving 
pseduopotentials method aided by B3LYP functional and an energy cutoff at 450 eV. 
The calculation is performed by using the spin polarized density functional theory*"^. 
For cell optimization, the supercell parameters are optimized in their charge states 



using PBE version of GGA with norm-conserving pseduopotentials. The geometry is 
further optimized by B3LYP method. The transition state search is implemented by 
the combination of linear synchronous transit and quadratic synchronous transit 
schemes*^ with 0.1 eV/A convergence threshold for the root mean square forces on 
the atoms^^'"^^'*^. K-point convergence test has been done and 5 irreducible k-points 
are used in sampling. For transition state search related calculations, only F point is 
used to reduce computational cost. Phonon dispersion and density of states cal- 
culation are computed by the finite displacement method implemented in CASTEP'**^. 
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